function plotrocks(filenumber)
%    N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
%time=readdata(0,6);
%     dtnom=readdata(0,7);
%dx=readdata(0,8);
%dz=readdata(0,9);
%mx=readdata(0,10);
%mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
%     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);
%     Nsurfaces=readdata(0,24);
%     surfaces=readdata(0,25);

load surface_time;
S=length(x);
timehistint=linspace(timehist(1),timehist(end),length(timehist)); dt=timehistint(2)-timehistint(1);

lowcut=0;
surfacehisthf=surfacehist;
for s=1:S
    %surfacehist=interp1(timehist,surfacehist,timehistint); surfacehist(1)=surfacehist(2);
    tempsurface=interp1(timehist,surfacehist(s,:),timehistint); tempsurface(1)=tempsurface(2);
    %temptemp=interp1(timehist,avtemphist(s,:),timehistint);
    %avtemphist(s,:)=bandpass(temptemp,lowcut,inf,dt);
    surfacehist(s,:)=bandpass(tempsurface,lowcut,inf,dt);
    %surfacehist(s,:)=surfacehist(s,:)-surfacehist(s,1);
    surfacehisthf(s,:)=bandpass(tempsurface,1/100,inf,dt);
    
    for d=1:length(depthrange)
        %temptemp=interp1(timehist,temphist(s,:,d),timehistint);
        %temptemp=temptemp-mean(temptemp);
        %temptemp=bandpass(temptemp,lowcut,inf,dt);
        %temphist(s,:,d)=temptemp;
    end
end
[TIMEHIST X]=meshgrid(timehistint,x);

for count=filenumber
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
    %     multigridlevels=readdata(count,4);
    %     bound=readdata(count,5);
    time=readdata(count,6);
    %     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    %     mcs=readdata(count,12);
    %     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
    %     vx=readdata(count,21);
    %     vz=readdata(count,22);
    %P=readdata(count,23);
    Nsurfaces=readdata(count,24);
    surfaces=readdata(count,25);

    x=zeros(1,S); z=zeros(1,T);
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end

    mTemp=readdata(count,15);
    [gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
    ind=find(gridtemp<1);
    for t=1:T
        gridtemp(:,t)=gridtemp(:,t)-mean(gridtemp(:,t));
    end
    gridtemp(ind)=nan;
    [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
    
    currentsurface=interp2(X',TIMEHIST',surfacehist',x,ones(1,S)*time/(60*60*24*365*1e6));
    
    figure(1); clf;
    subplot(2,2,[1 3]); pcolor(X,TIMEHIST,-surfacehist); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; hold on; plot([x(1) x(S)], [1 1]*time/(60*60*24*365*1e6),'k','linewidth',2);
    subplot(2,2,2); plot(x,-currentsurface); axis([0 2100e3 -170 170]); colorbar;
    subplot(2,2,4); pcolor(Xgrid,Zgrid,gridtemp); shading flat; set(gca,'ydir','rev'); axis equal; caxis([-150 150]); colorbar; axis([0 2000000 0 710000]); %title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    
    print(gcf, '-djpeg100','-r100',['..\outfiles\surf_surf_temp' num2str(count)]);

    
    
 
end
end